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INTRODUCTION 

The primary purpose of this paper is to describe a family of very 
powerful numerical integration techniques for use at the Lewis Research 
Center. The discussion and evaluation of these techniques will be of 
interest, however, at all computing centers where the techniques are not 
already receiving wide usage. Although the techniques are not new, or 
original, an appreciation of the savings in computor time that may be 
achieved by their application seems sadly lacking. In fact, a recent 
book‘d on "numerical methods and computor programming" discusses Simpson's 
Rule but not Gauss' formula. 

A comparison is made between these techniques and a modified 
Simpson's Rule program, STMPSlj' written by Drs. T. Fessler and W. Ford 
of the Lewis Research Center. A brief description of SIMPS1 and how it 
is related to ordinary Simpson's Rule will be presented herein. 

A description of the Fortran IV programs and their call vectors is 
given in the Appendix. The format presented herein is that of the Lewis 
Monitor Manual. These programs, with the exception of QUAD4, have been 
programmed by H. E. Renkel of the Lewis Research Center. 


GAUSSIAN QUADRATURE 


Gauss ' Formula 

According to Scarborough^, "The most accurate of the quadrature for- 
mulas in ordinary use is known as Gauss' formula." This formula can be 
written in the following form: 



N 


f (x)dx VW 


a-i. 


(i) 


The normalization to the region (-l,+l) imposes no real restriction (in 
fact, the normalization is effected internally in the FORTRAN programs). 
The abscissas x a and the weight coefficients have prescribed values 

for each N. Kopal^ discusses the evaluation of x and H for any 
given N. 


The reason Gauss's formula has not received more complete acceptance 
is given by Ralston 4 (1960): "Until recently Gaussian type quadrature 
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formulas were seldom used because in most cases the abscissas are 
irrational numbers, which, of course, makes them inconvenient for hand 
computation. However, on digital computers this argument has little 
weight . . . " 

Gauss' formula (eq. (l)) has been programmed in single and double 
precision. These function subroutines and their call vectors are de- 
scribed in the appendix. 


Weight Functions 


A series of formulas belonging to the Gaussian quadrature family is 
of the form 



W(x)f ' x) 


dx 


-£ 


- P; 


X. 


GO 


( 2 ) 


where W(x) is a given weight function. A list of the more common for- 
mulas of this type and their definition is given in table I. In the last 
column of table I the FORTRAN IV subprogram names are given. SQUAD1 and 
DSQUAD1 are subroutines. All the others are function subprograms. The 
prefix D denotes double precision programs. These subprograms are more 
fully described in the appendix. 


To appreciate the desirability of QUAD2 , it must be recognized that 
all commonly used integration formulas implicitly assume that the inte- 
grand may be approximated (over the region of integration) by a polynominal. 
Germane to the weight vs it -must also be recognized that no polynominal 
expansion in the independent variable can represent a curve that possesses 
an infinite derivative within the region of interest. 

'The comments for the case when the integrand has an infinite deriva- 
tive apply even more so when the integrand is singular in the region of 
integration (viz., QUAD 3, QUAD 4). 

Note that the formulas given herein are all open-ended formulas; that 
is, the integration formulas do net evaluate the integrand at the end- 
points. This property will be used later to advantage when the integrand 
has an implicit singularity. For further properties of Gaussian quadra- 
ture consult reference 5. 


The coefficients and abscissas of other weighted-Gauss ian quadrature 
formulas are no oe found in reference 5 and recent literature (such as the 
Journal of Mathematics of computation) . 


Multiple Integrations 

It is often necessary to numerically integrate equations of the fol- 
lowing type : 
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TABLE I. - GAUSSIAN QUADRATURE FORMULAS 



U 


Private communication from Peter Sockol. 
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dy dx 


(3) 


Subroutines SQUAD 1 and DSQUAD1 are presented in the appendix for this 
purpose. 


EVALUATION 


Usually, a reasonable indication of the efficiency of a method for 
numerical integration is given by the number of times the integrand must 
be evaluated to achieve a particular accuracy. Table II presents the re- 
sults of numerically integrating the integral. 

1 

W(x)f(x) dx (4) 


for various weights W(x) and functions f(x) by SIMPS1 and the appropri- 
ate Gaussian quadrature formula. The singular integrands (examples 3, 4, 

7 to 10, and 12 in table II'; cannot, of course, be evaluated by Simpson's 
Rule to the singularity • i here , zero) but can only be_evaluated to within 
some e of it. Results for e - 10“ 4 and e « 1.1 are presented in 
the table. Ng denotes the number of times SIMPS1 evaluated the integrand. 
The error, in examples 1 to 12, is given in units of the eighth digit, and 
represents the absolute error. The Gaussian quadrature errors are tabu- 
lated for 3, 5, 7, 9, and II integration points, respectively. The last 
column in the table shows the ratio of Ng to Nj where Nq represents 
the number of times the Gaussian quadrature needed to evaluate the inte- 
grand to achieve an error comparable to SIMES1. 

Examples 1 to 12 were chosen so that the integral could be performed 
exactly and an accurate absolute error could be obtained. Example 13 is 
an evaluation of the following integral: 



where 


. . in yl + a + -/l ± cpfjX dx 

yi + q^xT Va + 


Cp\ x} x f 5x" 
A - 0. 013 


( 5 ) 


This integral was encountered by the author in a recent problem and is 
included to show that even with a not too badly behaved integrand, the 
ratio Ng/Ng may exceed 10. 
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TABLE II. - RESULTS OF EVALUATING f W(x) f(x) dx 

JO 


Example 

W(x) 

f(x) 

SIMPS1 

Gaussian quadratures 

Ng 




% 

Error 

Routine 

Error 

% 







3 

5 

7 

9 

11 


1 

1 

1 + X 

5 

0 

QUAD1 

0 


0 

0 

-1 

1.6 

2 



49 

-10 

QUAD2 

0 


0 

0 

0 

16 

3 

W x 


a 117 

-4 

QUAD3 

-1 

0 

-1 

-1 

-1 

39 

4 

In x 


a 97 

1 

QUAD4 

0 

0 

1 

1 

1 

32 

5 

1 

1/(1 + x) 

17 

-1 

QUAD1 

- 

-4 

-1 

-2 

-2 

3.4 

6 

v* 


57 

-85 

QUAD2 

- 

-2 

-1 

-2 

-2 

17 

7 

i/V* 


a 137 

-2 

QUAD3 

- 

0 

0 

0 

0 

27 

8 

In x 


| a 109 

+8 

QUAD4 

- 

4 

5 

4 

6 

22 

9 

l/Vx 


^157 

-3 

QUAD3 

- 

0 

0 

0 

0 

31 

10 

In x 


b 117 

-3 

QUAD4 

- 

4 

5 

4 

6 

! 

23 

11 

1 

10e" 10x 

61 

-7 

QUAD1 

“ 

- 

-65 

-2 

-3 

7 

12 

l/n/x 

-Jx cos X 

b 169 

-25 

QUAD3 

31 

0 

0 

-2 


34 

■ 

-2 

13 

1 


98 

- 

QUAD1 

- 

- 

c 9 

C 0 

c o 

11 


(a) e = 10" 4 . 

(b) e = 10‘ 5 . 

(c) Error relative to SIMPS1. 
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Example 12 actually represents the integration of 



cos x 


2-y/sin 


dx 


( 6 ) 


where there is an implicit 1 /-^/x singularity in the integrand. This 
integral is easily evaluated by QUAD3 when expressed in the form 



cos x 1 




sin x 



dx 


where the function cos x/2-y/sin x/x is well behaved. Note that this 
integration does not necessitate a prior evaluation cf sin x/x at 
x = 0 since, as previously mentioned, the Gaussian quadratures do not 
evaluate the integrand at the end points. 

The range of values of the ratio N g /N G shown in table II must be 
considered from two different perspectives: as a measure of the efficiency 

of Gaussian quadrature methods compared with ordinary Simpson Integration, 
and as a measure of the Gaussian quadrature methods compared with SIMPS1 
(presumably only at lewis, to date). SIMPS1 is an ingenious function sub- 
program for the evaluation of integrals wherein the integrand may have many 
sharp peaks. SIMPS1 does not subdivide the region of integration uniformly, 
but preferentially subdivides only those subregions where it is necessary 
to achieve a prescribed maximum error. The basic integration technique 
employed by SIMPS1 is Simpson's Rule. Hence, SIMPS1, is a Simpson's Rule 
technique which minimizes the evaluation of the integrand. In other words, 
to obtain the same accuracy with a normal Simpson's Rule would require 
evaluation of the integrand at intervals equal to the smallest interval 
employed by SIMP31 over the whole range. This would in general, be many 
times the values shown for Ng in table II for evaluation of the same inte- 
grals. Therefore, the ratios Ng/Ng given in table II are very conservative 
estimates of the efficiency of Gaussian quadrature techniques when compared 
with ordinary Simpson's Rule. 

The ratio Ng/N. G does not provide a fair comparison of the Gaussian 
quadrature formulas with SIMPS1 since the latter may employ each evaluation 
of the integrand many times and use many more operations than the former. 
Table III shows a comparison of three examples from table II on the basis of 
time. The ratio Tg/'Tg shows the increase in speed available with the 
Gaussian quadrature, where Tg and Tg are the times needed to perform the 
integration by SIMPS1 and Gaussian quadrature, respectively. In the last 
column of table III, Tj represents the time spent in evaluating the inte- 
grand. Hence, the ratio represents a comparison of the time spent by the 
two methods in operating on the integrand evaluations. 


Table IV gives an indication of the accuracy obtainable with a double 
precision Gauss' formula (DQ'JADl). 


7 


TABLE III. - COMPARISON BETWEEN NUMBER OF EVALUATIONS 
OF INTEGRAND AND COMPUTATION TIME 


Example 

from 

table II 

, 

Evaluation of 
integrand 

Time 

Ts/Tq 

in 

EH 



n s /n g 

5 

17 

5 

3.4 

4.7 

5.3 

9 

157 

5 

31 

91 

78 

12 

169 

5 

22 

46 

177 


TABLE- IV. - DOUBLE PRECISION EVALUATION 
OF f (x)dx BY GAUSS' FORMULA 


f(x) 


N 


Absolute error 


<LO 


-15 


<LO' 


•15 
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CONCLUSIONS 

! It has "been shown that Gaussian quadrature is highly superior to 
ordinary Simpson's Rule or SIMPS1 for a rather wide class of functions. 
The examples considered show an increase in actual speed of the Gaussian 
quadrature routines over SIMPS1 of from '4. 7 to 91 times. Since the ma- 
jority of numerical integrations performed at large computing centers 
treat far less pathological integrands than those for which SIMPS1 was 
designed, it would appear advantageous to use the more rapid Gaussian 
quadratures whenever possible. One advantage of SIMPS1, however, is 
that it does provide an excellent means of determining the optimum 
number Nq needed for a given accuracy with Gaussian quadratures. 
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APPENDIX - FORTRAN PROGRAM DESCRIPTIONS 
QUADl/lNTEGRATION 

Purpose Function subprogram used to perform the numerical integration of 

~ — nXF 

/ f(x) dx 
JX 0 

"by Gauss' formula. A very powerful method for minimizing machine 
time in production runs when the integrand and first derivative 
are continuous in the region of integration. 

Calling 

Sequence QUAD1(N,NS,X0,XF,F0FX:) 

N is the number of points per division. 

NS is the number of divisions into which XF-XO is divided (this 
is provided in case 16 points is not sufficient - see below) 

FOFX is an external function subprogram. 

N = 3(1)16* 9 D 


Source 

Language FORTRAN TV 

Reference 8 


IQUAD1 /DOUBLE PRECISION INTEGRATION 
Purpose Double Precision version of QUAD1 (q. v. ) 


Calling 

Sequence DQUAD1(N,NS,XC,XF,F0FX) 

FOFX must be a double precision function 
N = 3(1)15, 15 D 

N = 16(4)24(8)48(16)96, 17 D** 

Reference 8 


-X- 

This notation denotes that N takes on all values from 3 to 16 at 
unit increments; for example, N = 3(1)6(2)12 signifies N = 3,4,5,6,8,10,12. 

Data available to 21 D on punched cards. 
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Purpose 


Calling 

Sequence 


EXAMPLE 


Purpose 


SQUADl/lHTEGRATION OF MULTIPLE INTEGRALS 


Subroutine subprogram used to perform the numerical integra- 
tions of the types 



g(x,z) 


dz dx 


by repeated application of Gauss' formula. (cf. QUADl) 


SQUAD1 (MODE,N,XO,XF,X,Y,ANSWR) 

MODE = 1, Obtain arguments X(N). 

= 2, Use evaluations of integrand Y(X(N)) to obtain ANSWR. 
N number of points of integration 
X arguments - must be dimensioned array name X(N) 

Y - must be dimensioned array name Y(N) 

XNSWR result of integration 



G(x') dx' dx 


A,B, C constants 

F(X), G(X') function subprograms 


Dimension XF(NF) , YF(NF) ,XG(NG) , YG(NG) 


CALL SQUAD1 (1,NF,B,C,XF, YF,FINTGL) 

DO 10 I = 1,NF 

CALL SQUAD1 (l,NG,A,XF(l) ,XG, YG,ANSWR) 
DO 20 J = 1,NG 
20 YG( J) = G(XG( J) ) / 

CALL SQUAD1 (2,NG,A,XF(l),XG,YG,AMSWR) 

10 YF(l) =? ARSWR* F( XF( I ) ) 

CALL SQUAD1 ( 2,MF,B, C,XF, YF,FINTGL) 


DSQUAD1 /DOUBLE PRECISION MULTIPLE INTEGRATION 
Double precision version of SQUAD1 (q.v. ) 
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QUAD2 /INTEGRATION 


Purpose 


Calling 

Sequence 


Function subprogram used to perform the numerical integration of 



-/ x - XO f (x) dx 


by weighted Gaussian quadrature. A very powerful technique for 
minimizing the machine time in production runs when f(x) and 
its first derivative are continuous throughout the region of 
integration. 


QUAD2(N,X0,XF,F0FX) 
See QUAD1 
N = 2(1)8 8 D 


Reference 5 


Purpose 


QUAD3/INTEGRATI0N 


Function subprogram. Used to perform the numerical integration 
of 


P XF . 

/ „ f(x) 

Jx o 7F - X0 1 


dx 


by weighted Gaussian quadrature. A very powerful method for 
minimizing the machine time in production runs when f(x) and 
its first derivate are continuous throughout the region of 
integration. 


Calling 

Sequence QUAD3(N,X0,XF,F0FX) 

See QUAD1 

N = 2(1)8(2)12(4)24(8)48 17 D 

Reference 6 


DQUAD5/D0UBLE PRECISION INTEGRATION 
Purpose Double precision version of QUAD3 (q.v. ) 
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Purpose 


Calling 

Sequence 


Method 


'Other 

Subroutines 

Employed 

Reference 


QUAD4/ IMEGRATI ON 

Function subprogram. Used to perform numerical integration 
of 


f(x) In | x - X0| dx 

by weighted Gaussian quadrature. A very powerful method for 
minimizing the machine time in production runs when f(x) and 
its first derivative are continuous throughout the region of 
integration . 



QUAD4( N, XD, XF, FOFX ) 
(See QUAD1 
N = 3(1)11 9 D 



f ( x ) In 


- X0| dx = AX 



f(AC • y + XO) In y dy 


+ In 



f(x) dx 


where 


AX = XF - XO 


y = 


x - XO 


AX 


QUAD1 (q.v.) 

Private Communication from Peter Sockol 
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QUAD5 /INTEGRATION 

Purpose Function subprogram used to perform the numerical integration 

of 


e“ x f(x) dx 


by the Gauss -Laguerre quadrature formula. A very powerful 
method for minimizing the machine time in production runs when 
f(x) and its first derivative are continuous throughout the 
region of integration. 



Calling 

Sequence QUAD5 ( N, FOFX ) 

See QUAD1 

Reference 7 


Purpose 


QUAD6 /INTEGRATION 

Function subprogram used to perform the numerical integration 
of 


n +0 ° 2 

r e~ x f(x) dx 

La- oo 

by the Gauss -Hermite quadrature formula. A very powerful 
method for minimizing machine times in production runs when 
f(x) and its first derivative are continuous throughout the 
range of integration. 


Calling 

Sequence QUAD6(N,F0FX) 
See QJJAD1 


Reference 


9 
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